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SECTION  1 


INTRODUCTION 

Panel  flutter  is  a  self-induced  oscillation  of  a  thin- 
walled  structure  caused  by  an  increase  of  aerodynamic  forces 
resulting  from  panel  deformation  in  an  airstream.  Panel  flutter 
does  differ  from  wing  flutter  in  that  the  nature  of  panel  flutter 
is  generally  not  as  catastrophic  as  wing  flutter.  Linear  flutter 
prediction  techniques  allow  one  to  determine  the  locus  of 
neutrally  stable  oscillations  which  is  called  the  flutter  boun¬ 
dary.  Nonlinear  panel  flutter;  however,  is  characterized  by  a 
periodic  oscillation  of  finite  amplitude  at  and  above  the  linear 
flutter  boundary.  Presently,  the  two  methods  available  for  non¬ 
linear  flutter  prediction  are  theoretical  analysis  and  experi¬ 
mental  determination. 

At  the  present  time  there  is  poor  correlation  between 
experimental  and  theoretical  flutter  prediction  in  the  transonic 
Mach  number  regime.  In  the  past,  the  theoretical  prediction  of 
panel  flutter  has  been  overly  conservative  due  to  simplifying 
assumptions  made  because  of  the  complexity  of  panel  flutter. 
However,  these  assumptions  are  usually  so  restrictive  that  there 
is  great  variance  between  experimental  results  and  theoretical 
predictions  at  transonic  speeds. 

The  theoretical  studies  have  been  restricted  by  use  of  a 
simplified  linear  aerodynamic  theory  and  an  imprecise  idealiza¬ 
tion  of  the  panel  support  conditions  which  result  when  the 
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nonlinear  midplane  stresses  (i.e.,  use  of  small  amplitude 
oscillation)  are  neglected.  This  study  removes  the  above  two 
simplifying  assumptions  in  that  the  large  deflection  equations  of 
a  panel  in  a  nonlinear  transonic  airstream  with  shocks  present, 
are  formulated  for  solution. 

Dowell*  has  presented  a  study  of  the  nonlinear  flutter  of  a 
flat  panel  using  a  Galerkin  procedure  wherein  he  numerically 
integrates  the  resulting  ordinary  differential  equation  of  motion 
in  the  time  variable.  Dowell  has  also  used  a  nonlinear  plate 
theory  but  a  linear  aerodynamic  theory  in  the  form  of  both  a 
quasi-steady  and  a  full  unsteady  theory.  The  availability  of 
such  an  analysis  permits  the  consideration  of  the  panel  post¬ 
flutter  behavior  only  over  the  subsonic  and  supersonic  Mach 
number  range.  In  the  transonic  range  it  is  necessary  to  use  a 
nonlinear  aerodynamic  theory  since  the  governing  differential 
equation  is  nonlinear  and  contains  both  subsonic  and  supersonic 
regions  together  with  shocks. 

The  method  to  be  used  herein  will  be  briefly  outlined.  The 
equations  of  motion  for  the  transverse  oscillation  of  a  panel  are 
obtained.  It  is  assumed  that  the  large  deformations  can  be  ade¬ 
quately  described  by  use  of  the  Von-Karman  large  strain-displace¬ 
ment  relations.  The  flutter  motion  of  the  panel  is  described  by 
a  simple  two  mode  approximation.  Galerkin's  method  is  used  to 
obtain  the  differential  equation  in  time  governing  the  panel 
response  to  a  time  dependent  aerodynamic  pressure  field.  The 
aerodynamic  pressure  is  to  be  transonic  and  as  such  is  governed 
by  the  small-disturbance,  potential  transonic  flow  differential 
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equation.  The  solution  to  the  nonlinear  transonic  differential 
equation  is  obtained  using  a  finite  difference  method  developed 
by  Ballhaus  and  Goorjian^  and  modified  for  large  values  of 
reduced  frequency  by  Huizing  and  van  der  Vooren^.  This  time 
dependent  response  investigation  requires  simultaneous  integra¬ 
tion  of  the  panel  large  deflection  equation  and  the  transonic 
aerodynamic  equations  in  time. 

Objectives 

The  main  objective  of  this  investigation  is  to  determine  a 
procedure  for  the  prediction  of  the  large  deflection  response  of 
a  structural  panel  situated  in  a  transonic  airstream.  Von-Karman 
large  deflection  equations  are  to  be  integrated  simultaneously 
with  solutions  of  the  nonlinear  transonic  equations  obtained  from 
a  finite  difference,  alternating  direction  implicit  scheme 
referred  to  as  LTRAN2 .  The  prime  task  then  is  to  incorporate 
this  integration  scheme  into  a  modified  version  of  LTRAN2  called 
NLR-LTRAN2.  This  procedure  will  yield  theoretical  predictions  of 
flutter  speed  which  should  compare  with  experimentally  determined 
speeds  within  the  transonic  Mach  number  region. 

Large  Deflection  Equations  for  Panel  Oscillation 

Consider  the  isotropic  thin-walled  panel  shown  in  Figure  1. 
It  is  exposed  to  a  transonic  airstream  in  which  shocks  may  occur 
and  is  subjected  to  a  longitudinal  axial  load.  It  is  assumed 
that  the  large  oscillations  of  the  panel  in  the  presence  of  mid¬ 
surface  stresses  are  governed  by  the  nonlinear  Von-Karman  strain 
displacement  relations.  These  relations  are  given  as: 
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E  =  10  x  106  lb/in2 
h  =  0.03  in 
a  =  b  =  10  in 


Plate  edges  are  simply 


sorted 


Figure  1.  A  Thin  Uniform  Elastic  Plate  Exposed 
to  a  Transonic  Airstream 
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where  z  is  measured  from  the  midsurface  of  the  panel.  Equation 
(1)  is  a  modification  of  a  small  deflection  theory  to  include  the 
first  order  effects  of  midsurface  stretching  necessary  to  inves¬ 
tigate  large  deflections. 

Now  Hamilton's  variational  principle  is  enforced  to  obtain 
the  partial  differential  equations  of  panel  motion.  For  a  three- 
dimensional  plate,  Von-Karman's  large  deflection  equations  (see 
Eastep  and  McIntosh5  for  a  derivation)  are: 


D,Vh  r? 

3y  3x 


2  2  2  2  2  2 
3  F  3  w  3  F  3  w  -  3  w  3  w 

3x3y  3x3y  +2  ,2x2,2 

1  1  3x  3y  3  x  3  x 


2  2 

,  M  3  w  U  3  w  ,  .  ,  .  . 

+  N  =■  =  -ph  =■  +  A(w,t) 

y  ai.  ^ 


3y 


3 1 


and 


(2) 


(3) 


where  w  is  the  plate  transverse  deflection  while  F  is  the  Airy 
stress  function.  The  aerodynamic  pressure  loading,  A,  is  the 
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increased  incremental  aerodynamic  pressure  caused  by  the  panel 
deformation  and  is  given  by: 


A(w,t)  =  -P„[j£  +  UQ  f£] 


where  the  velocity  potential  4  must  satisfy  the  small  disturbance 
transonic  differential  equation: 


(lV)  ^4  -  (Y  ♦  1)  ll  -  2*1  -  2M-2 


- '  2 
3x 


+  1^|  =  0 
ay 


3*  3x2  -  3t2 


3t3x 


The  numerical  computation  method  used  for  determining  the 
unsteady  transonic  flow  field  is  based  on  the  alternating  direc¬ 
tion  implicit  (ADI)3  procedure  modified  to  consider  flows  of 
moderately  high  reduced  frequencies  and  panel  deformations.  The 
procedure  is  referred  to  as  NLR-LTRAN2.4 

The  system  of  Equations  (2)  through  (5)  will  be  solved  by 
the  Galerkin  method  for  a  simply  supported  plate  undergoing 
cylindrical  bending.  That  is: 


W(x,t)  =  h[c.(t)sin  7^  +  C_(t)sin  —7^] 

X  a  4  a 


I-  I 


SECTION  2 


APPROACH 


Two-Dimensional  Flow 


Consider  for  simplification  the  flow  about  a  two- 
dimensional  wing  panel  of  given  thickness  undergoing  an  assumed 
chordwise  deformation  and  the  prediction  of  the  oscillatory  aero¬ 
dynamic  pressures  in  the  transonic  speed  regime.  Figure  2  shows 
a  panel  that  is  simple  supported  at  both  ends.  For  this  study 
the  aerodynamic  loading  was  assumed  to  act  only  on  one  side  of 
the  panel.  The  undeformed  shape  of  the  panel  was  represented  by 
a  half-sine  wave  for  both  the  upper  and  lower  surfaces.  With  a 
chord  of  10  inches,  the  thickness  to  chord  ratio  of  the  panel  was 
0.003. 

For  two-dimensional  flow,  all  derivatives  with  respect  to  y 
are  zero.  Therefore,  th  Von-Karman  large  deflection  plate  equa¬ 
tion  becomes: 


4  2  2 

^  3  w  -  3  w  .  y,  3  w 

aw4  x  3x2  m  3 1 2 


=  -  (P  -  P-) 


where  the  longitudinal  axial  load  Nx  is  defined  as: 


■,-!5  l  ‘H’2** 


Galerkin's  Method 


The  procedure  used  in  solving  Equation  (7)  was  to  apply 
Galerkin's  Method  and  a  step-by-step  time  integration  method. 
Galerkin's  Method  is  an  assumed  mode  technique  which  reduces  a 
partial  differential  equation  with  independent  variables  x  and  t 
to  a  set  of  simultaneous  ordinary  differential  equations  with 
independent  variable  t. 

For  the  Galerkin  procedure,  a  series  of  displacment  func¬ 
tions  which  satisfy  both  the  geometric  and  force  boundary  condi¬ 
tions  are  assumed.  In  practice,  the  more  terms  used  in  the 
series,  the  more  accurate  the  answer  will  become.  For  the  simple 
supported  beam,  the  geometric  boundary  conditions  imposed  on  the 
problem  are: 

w(x,t)  |  ®  *  0,  for  all  time  (8) 


and  the  force  boundary  conditions  are: 


El 


3  w( x, t ) 
3x2 


a 

o 


0,  for  all  time 


(9) 


The  transverse  deflection  w(x,t)  can  then  be  approximated 
by  a  series  of  displacement  functions  weighted  by  unknown  coeffi¬ 
cients  of  time.  Therefore, 


N 

w(x,t)  *  w( x , t )  *  h  l  C.(t)4>.(x),  (10a) 

i-1  1  1 


or 
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w(x,t)  =  h(C ,  ( t )  sin  ^  +  C-(t)sin  + 

1  a  4  d 


,  „  , .  .  Nirx-, 
+  CN  t  > sm  — ) 


where  the  displacement  functions  are  defined  as: 


,  .  lirx 

<|> .  (x)  =  sin  — — 

1  a 


(10b) 


(10c) 


The  basic  procedure  in  applying  the  Galerkin  Method  is  to 
substitute  Equation  (10a)  into  Equation  (7);  multiply  the  resul¬ 
ting  equation  by  <j>^(x),  i  =  1,...N;  and  integrate  over  the  domain 
of  the  problem.  This  can  be  expressed  in  equation  form  as: 


a  4-  2-  2— 

f  fun  3  w  .  r.  3  w  .2  3  w 

/  (hD  — y  -  hN  — 2  +  pmh  — 2  + 
o  3X*  X  3w^  m  at* 


h(P  -  P«))  sin  — —  dx  =  0  for  i  ■  1,...N 

di 


This  will  result  in  N  simultaneous  ordinary  differential  equa¬ 
tions  with  unknowns  C^(t). 

Piston  Theory  Aerodynamics 

To  become  familiar  with  using  the  Von-Karman  large  dis¬ 
placement  equations  and  to  check  the  resulting  equations  after 
applying  Galerkin's  Method,  a  simple  aerodynamic  theory  was  used 
as  the  forcing  function.  For  this  study,  Piston  Theory  aero¬ 
dynamics  was  used  because  of  the  ease  in  which  (P  -  P<»)  can  be 
calculated,  and  secondly  because  of  the  availability  of  previous 
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•  1  2 
studies  conducted  by  Prof.  Dowell  '.  In  the  development  of 

Piston  Theory,  the  flow  was  assumed  to  be  highly  supersonic  or: 

0  =  /  M2  -  1  =  M  (12) 

CD  OD 

With  this  assumption,  the  pressure  distribution  could  be 
expressed  as: 

,p  -  p->  -  |a  [|f  +  i|fj  (i3) 

Time  Integration  Procedure 

After  applying  Galerkin's  Method  and  transforming  the  equa¬ 
tion  into  a  nondimensional  form,  the  set  of  simultaneous  differ¬ 
ential  equations  can  be  obtained  from  Equation  (11): 


[Ml  {U}  +  [C]  {U}  +  [K]  { U}  =  {P} 


(14) 


where  [M] ,  [C] ,  and  [K]  are  the  mass,  damping  and  stiffness 
matrices  respectively,  {U}  is  the  vector  of  the  unknown  coeffi¬ 
cients  of  time,  and  {P>  is  a  vector  of  the  aerodynamic  loads. 

The  time  derivatives  are  with  respect  to  the  nondimensional  time. 

The  nondimensional  variables  used  in  obtaining  Equation 
(14)  include: 


X 


2qa3 

~w  7  v 


au 

V 


T 


Ut 


(15) 
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The  solution  of  Equation  (14)  was  obtained  using  a  step-by- 
step  time  integration  finite  difference  approach6.  Assuming  a 
linear  variation  of  acceleration,  the  velocities  and  displace¬ 
ments  at  the  end  of  a  small  time  interval  can  be  expressed  as: 


(U), 

(U>t 


+T-  + 


(16) 

(17) 


where  At  is  the  time  step  and  t-At  is  the  previous  time.  Substi¬ 
tution  of  Equations  (16)  and  (17)  into  Equation  (14)  gives  the 
acceleration  at  the  new  time  step  as: 


(U>T  *  CFI  [  CP}  t  -  (CJ  {v}  -  CK]{w}] 


(18) 


where 


(F]  =  [  (M]  +  [C]  +  [K]  J 


(19) 


{V>  "  {6>t-At  +-T  ^>t-At 


(20) 


(w>  =  (U)t_At  -  A*{U}t_At  +  {“U}t-Ar 


(21) 


Equation  (18)  can  then  be  used  to  find  the  velocity  (Equation  16) 
and  displacement  Equation  (17)  at  the  new  time  step. 


The  vector  {P}  is  obtained  numerically  by  solving  the 
governing  aerodynamic  equations.  The  displacment  and  velocities 


needed  for  computing  {P}  are  based  on  the  values  obtained  at  the 
old  time  tep.  The  time  step  was  chosen  small  enough  such  that  no 
numerical  instabilities  would  occur  during  the  solution.  In 
general,  {P}  can  be  shown  to  be: 


{P}  = 


-(£>"* 


K2h 


7  (2=2=)  sin  ,  i  Si 
0  q  a  a 
r 


/  [— )  "in  2"  7  —  } 
}  (EjEi)  sln  „,1SS 

0  q  a  a 


(22) 


where  (P  -  P«)  can  be  determind  from  either  Piston  Theory  or 
LTRAN2. 

Transonic  Aerodynamic  Pressures 
The  problem  of  interest  is  the  flow  about  a  wing  panel  of 
given  thickness  undergoing  an  assumed  chordwise  deformation 
(herein  selected  to  be  two  modes)  and  the  prediction  of  the 
oscillatory  aerodynamic  pressures  in  the  transonic  speed  regime. 
The  numerical  computation  method  used  for  the  unsteady  flow  field 

is  based  on  the  alternating  direction  implicit  (ADI)  procedure 

*1 

first  formulated  by  Ballhaus  and  Goorjian  .  The  procedure  was 
modified  by  van  der  Vooren^  to  consider  flows  of  moderately  high 
reduced  frequencies  and  panel  deformations.  The  small  distur¬ 
bance  velocity  potential  equation  can  be  derived  for  transonic 
flow  by  assuming  an  inviscid  isentropic  fluid  with  only  weak 
shocks  existing.  The  resulting  moderately  high  reduced  frequency 
differential  equation  solved  in  the  code  NLR-LTRAN2  is  Equation 
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>•  Ktt* 


(5)  with  the  $  term  dropped.  In  addition  NLR-LTRAN2  retained 
necessary  unsteady  terms  in  the  boundary  conditions  on  the  panel 
and  the  wake  condition  of  zero  pressure  differential  in  the 
wake.  These  modifications  of  the  original  LTRAN2  code  of 
Reference  3  allowed  reduced  frequencies  up  to  k  =  0.8  to  be  con¬ 
sidered.  The  unsteady  two-dimensional,  transonic  small  distur¬ 
bance  equation  solved  in  the  code  NLR-LTRAN2  is: 

-2kM  2<(.  +  [  1— M  2  -  (y+1)M  2 4>  U  +  *  =0  (23) 

where  <j>(x,y,t)  is  the  small  disturbance  velocity  potential 
resulting  from  the  panel  deformation  w(x,t).  The  boundary  con¬ 
ditions  which  must  be  satisfied  are: 

<k  (x,y=0,t)  =  "r~  +  —r-  on  the  panel  surface  (24) 

y  d  X  d  L 

and 

ACp(x,y=0,t)  =  A(<|>x  +  =  0  across  the  wake  (25) 

where  w(x,t)  is  the  panel  deformation.  At  large  distances  away 
from  the  panel  we  require  that: 

♦x2  +  *y2  ♦  0  (26) 

For  a  prescribed  panel  deformation  w(x,t)  then  NLR-LTRAN2  can  be 
used  to  determine  the  velocity  potential  from  Equation  (23)  using 
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the  AID  scheme  of  Reference  3.  With  the  velocity  potential  thus 
determined/  the  pressure  field  and  hence  the  pressure  difference 
of  the  panel  can  be  obtained  from: 

Cp  =  ‘2Ux  +  Mt>  (27) 

The  panel  considered  in  this  study  is  a  symmetric  section 
composed  of  sinusoidal  arcs  of  max  thickness  ratio  as  shown  in 
Figure  2.  This  panel  is  placed  at  zero  angle  of  incidence  in  a 
transonic  airstream  of  Mach  number  0.85.  The  chord  deformation 
is  selected  by  the  two-mode  assumption  as  given  in  Equation  (10) 
where  the  maximum  slope  of  the  chord  deformation  is  limited  to  1 
degree.  The  code  NLR-TRAN2  is  used  to  calculate  the  steady  state 
initial  conditions  for  the  10%  thick  panel  and  plot  of  the  upper 
surface  coefficient  of  pressure  shown  in  Figure  3.  The 
occurrence  of  the  shock  wave  can  be  detected  at  the  70%  chord 
location.  In  addition,  shown  in  Figure  4,  is  the  upper  surface 
coefficient  of  pressure  for  the  first  mode  at  nondimensional  time 
of  4.  The  shock  wave  has  now  moved  to  the  new  location  of  90% 
chord  location.  For  assumed  values  of  the  modal  coefficients  C-^ 
and  C2  then  the  transonic  aerodynamic  pressure,  P,  required  in 
Equation  (22)  can  be  obtained  from  NLR-LTRAN2.  However,  in 
general,  the  determination  of  the  modal  coefficients  must  result 
from  a  simultaneous  integration  of  the  panel  equations  and  the 
transonic  aerodynamic  equations  as  described  in  the  previous 
section  on  the  time  integration  procedure. 
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SECTION  3 


RESULTS 


The  results  are  presented  as  two  tasks;  the  first  task 
involves  the  use  of  the  Von-Karman  large  deflection  equations  and 
Piston  Theory  aerodynamics.  The  second  task  involves  the  use  of 
two  single-degree-of-f reedom  models  being  loaded  by  transonic 
aerodynamics . 

Large  Displacement  Equations  and  Piston  Theory 
Assuming  a  two  mode  series  solution  for  the  Galerkin  proce¬ 
dure  (with  no  damping,  [C]  =  0),  the  mass  and  stiffness  matrices 
of  Equation  (19)  are: 

.  fi  ol 


(Ml  = 


(28) 


0  1 


and 


[K] 


k2xe 


3( Cj 2  +  4C22)  0 

0  16  +  22(C12  +  4C22) 


(29) 


The  mass  became  the  identity  matrix  as  a  result  of  the  nondimen- 
sional  form  of  the  final  equations  of  motion.  Also,  the  stiff¬ 
ness  is  nonlinear  since  it  is  a  function  of  the  square  of  the 
panel  displacements  and  C2. 

The  value  (P  -  P«)  can  be  found  as  a  function  of  the  panel 
displacement  and  velocities.  The  variation  of  (P  -  P~)  across 


the  chord  of  the  panel  can  then  be  weighted  by  each  of  the 
assumed  modes  and  integrated  to  obtain  the  force  vector. 

Assuming  a  two-term  representation  for  the  panel  displace¬ 
ment  and  using  Piston  Theory,  the  panel  deflection  at  limit  cycle 
for  increasing  dynamic  pressure  is  shown  in  Figure  5  as  a  solid 
line.  The  dash  lines  represent  similar  calculations  performed  by 
Prof.  Dowell  at  Princeton  several  years  earlier.  Point  solutions 
were  also  made  with  a  four-term  and  six-term  representation  of 
the  panel  displacement.  For  all  cases,  the  present  analysis 
correlates  very  well  with  previous  calculations. 

The  panel  mode  shape  is  a  function  of  the  aerodynamic  load¬ 
ing  and  is  presented  in  Figure  6.  Once  again,  the  data  agrees 
very  well  with  Dowell's  results.  The  maximum  deflection  of  the 
panel  was  shown  to  occur  near  the  70%  chord  point  for  all  varia¬ 
tions  in  the  aerodynamic  loading. 

Aeroelastic  Response  of  a  Single  Degree  of  Freedom  Airfoil 

For  this  task,  the  aeroelastic  response  of  an  airfoil  con¬ 
strained  to  deflect  in  either  pitch  or  plunge  was  determined 
using  the  NLR-LTRAN24  computer  code.  The  equation  of  motion  for 
an  airfoil  pitching  about  its  midchord  is: 

8C 

a  +  A^a  +  j  (30) 

ttu'K 

d 

For  this  case,  A^  =  .5,  A2  *  1.5  and,  p '  =  1000.  These  data 
were  selected  so  that  the  predictions  could  be  correlated  with 
similar  results  found  in  Reference  7. 
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Figure  6.  Panel  Mode  Shape  at  Limit  Cycle 
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The  dynamic  response  o£  the  flat  plate  pitching  about  the 
mid-chord  at  a  Mach  number  of  0.7  was  obtained  for  a  reduced  fre¬ 
quency  of  0.1.  The  aerodynamic  equation  was  integrated  in  time 
for  two  cycles  by  forcing  a  sinusoidal  variation  of  pitching 
angle  with  amplitude  of  0.01  radians.  The  free  motion  was 
started  at  the  end  of  the  second  cycle.  The  pitching  moment 
coefficient,  Cm,  was  determined  by  LTRAN2  at  each  time  interval 
and  used  in  Equation  (30)  to  determine  the  time  history.  A  con¬ 
verging  type  response  was  obtained  for  both  the  pitching  angle 
and  the  pitching  moment  indicating  that  the  flight  conditions  are 
below  the  instability  speed.  The  results  provided  in  Figure  7 
compared  very  well  with  data  from  Reference  7. 

Similarly,  for  a  flat  plate  plunging,  the  equation  of 
motion  is: 

-2C 

6  +  B^6  +  ^  (31) 

it  U  '  K 

a 

To  obtain  data  for  correlation  with  Reference  7,  the  coefficients 
B^,  B2  and  y '  were  selected  to  be  0.,  1.0  and  100  respec¬ 
tively.  Under  similar  starting  conditions  as  discussed  for  the 
pitching  airfoil  case,  the  plunging  response  and  the  coefficient 
of  lift  variation  with  time  is  presented  in  Figure  8  for  a  Mach 
number  of  0.7  and  a  reduced  frequency  of  0.1.  Both  the  plunging 
displacement  and  lift  coefficient  agreed  very  well  with  the 
Reference  7  results. 


SECTION  4 


CONCLUSIONS  AND  RECOMMENDATIONS 
Flutter  analyses  have  been  conducted  on  a  simply  supported 
panel  to  demonstrate  the  successful  combining  of  panel  Von  Karman 
large  deflection  equations  with  a  simple  linear  aerodynamic 
theory  (Piston)  for  determining  panel  response.  The  panel  re¬ 
sponse  was  determined  from  a  numerical  time  integration  scheme 
which  reproduced  the  results  presented  previously  by  Dowell.  In 
addition,  the  time  integration  scheme  was  successfully  used  to 
insure  the  simultaneous  integration  of  a  set  of  linear  structural 
equations  and  nonlinear  aerodynamic  equations.  Here,  the  simul¬ 
taneous  integration  scheme  was  used  to  determine  the  aeroelastic 
response  of  the  linear  pitching  or  plunging  of  an  airfoil  to  the 
nonlinear  aerodynamic  lift  and  moment  obtained  from  the  NLR- 
LTRAN2  computer  code.  The  response  obtained  for  pitching  or 
plunging  compared  favorably  to  those  obtained  by  Yang,  et  al.^ 

An  analysis  was  also  initiated  to  study  the  panel  response 
represented  by  a  nonlinear  structural  equation  (Von  Karman's 
large  deflection)  and  a  nonlinear  aerodynamic  equation  repre¬ 
sented  by  NLR-LTAN2 .  Here,  the  panel  response  was  represented 
with  an  assumed  modal  series  combination  even  though  both  the 
structural  equations  and  aerodynamic  equations  are  nonlinear. 
However,  to  date,  numerical  instabilities  with  the  computer  code 
NLR-LTRAN2  have  hampered  the  progress  and  obtaining  of  final 
results . 
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It  is  recommended  that  the  analysis  of  the  nonlinear 
response  of  a  panel  in  a  transonic  airstream  be  continued  using 
the  simultaneous  integration  of  structural  and  aerodynamic 
equations  described  herein.  Further,  it  is  recommended  that  the 
assumed  modal  technique  used  here  with  nonlinear  equations  be 
investigated.  Toward  this  end,  the  replacement  of  the  Von  Karman 
large  deflection  equations  with  a  finite-element  representation 
is  described  in  Appendix  A.  Unfortunately,  only  Piston  Theory 
aerodynamics  have  been  investigated  with  the  finite  element  model 
so  it  suggested  that  the  nonlinear  aerodynamic  theory  obtained 
from  NLR— LTRAN2  be  coupled  with  the  element  model  of  Appendix  A 
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APPENDIX  A 


PANEL  RESPONSE  OF  A  FINITE  ELEMENT  MODEL 
AT  POST  FLUTTER  CONDITIONS  USING  PISTON  THEORY 

Introduction 

Panel  flutter  is  a  self-induced  oscillation  of  a  thin- 
walled  structure  caused  by  increasing  aerodynamic  pressures 
resulting  from  the  panel  deformation  in  an  airstream.  Nonlinear 
panel  flutter  is  characterized  by  a  periodic  oscillation  of 
finite  amplitude  commonly  referred  to  as  limit  cycle.  Theore¬ 
tical  studies  in  the  past  have  been  restricted  by  the  use  of  a 
simplified  linear  aerodynamic  theory  and  an  inaccurate  repre¬ 
sentation  of  the  panel  support  conditions  which  result  when  the 
nonlinear  midplane  stresses  are  neglected  (assumption  of  small 
amplitudes).  References  1,  2  and  5  present  some  of  the  large 
amount  of  work  completed  in  this  technical  area.  Some  details 
have  been  extracted  from  these  references  for  use  in  the  study 
presented  here. 

The  main  objective  of  this  investigation  is  to  determine  a 
procedure  for  the  prediction  of  the  large  deflection  response  of 
structural  panels  in  a  transonic  airstream.  The  Von  Kartnan  large 
deflection  equations  are  to  be  integrated  simultaneously  with  the 
solutions  of  the  nonlinear  transonic  equations.  Since  this 
approach  requires  an  extremely  large  amount  of  computation,  the 
Von  Karman  large  deflection  equations  will  be  used  in  conjunction 
with  Piston  Theory  aerodynamics.  This  section  presents  only  the 
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results  of  the  calculations  with  Piston  Theory  and  a  finite 
element  model. 

The  thin  two-dimensional  plate  shown  in  Figure  A-l  has  been 
selected  as  the  configuration  to  be  investigated  in  this  study. 
The  transverse  aerodynamic  load  acting  on  only  one  side  of  the 
panel  is  of  known  spatial  and  temporal  distribution.  A  finite 
element  representation  was  selected  to  eliminate  the  modal 
approximations  necessary  in  the  previous  investigation  because  of 
the  nonlinear  Von  Karman  large  deflection  equation.  The  plate  is 
represented  as  consisting  of  twin  variable  node  plane-strain 
elements  each  of  unit  length.  A  thickness  of  0.03  inches  and  the 
material  properties  of  aluminum  are  used  because  of  studies  com¬ 
pleted  on  similar  plates  in  References  1,  2  and  5.  The  plate 
consists  of  54  nodes;  nodes  #1  and  #44  are  simply  supported  (con¬ 
strained  in  the  _x_»  X  and  X  directions),  while  the  other  nodes  are 
only  constrained  from  moving  in  the  _y_  direction.  Node  #54  repre¬ 
sents  the  three-quarter  span  of  the  plate. 

It  is  assumed  that  the  large  displacements  of  the  plate  are 
governed  by  the  nonlinear  Von  Karman  large  deflection  equations. 
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and  from  the  compatibility  equation: 


where  _w_  is  the  plate  transverse  deflection  and  _F_  is  the  Airy 
stress  function.  The  aerodynamic  pressure  loading,  A(x,t),  is 
the  increased  incremental  aerodynamic  pressure  caused  by  the 
plate  displacement.  As  mentioned  previously,  A(x,t)  for  this 
study  is  obtained  from  Piston  Theory.  These  equations  are  solved 
in  a  computer  code  referred  to  as  MAGNA,  Materially  Meometr^- 
cally  _Npnlinear  Ana^ys:*-s^ '  ^  using  a  finite  element  procedure. 

Discussion 

This  section  briefly  describes  the  MAGNA  computer  program 
capabilities,  and  presents  a  derivation  of  a  subroutine  for  MAGNA 
for  calculating  Piston  Theory  aerodynamic  pressures. 

MAGNA  Computer  Program 

The  MAGNA  computer  program  is  a  large  scale,  general  pur¬ 
pose  finite  element  system  intended  for  the  nonlinear  (large 
deflection)  analysis  of  complex  engineering  structures.  MAGNA 
has  been  developed  primarily  for  the  efficient  solution  of  three- 
dimensional  problems  involving  many  degrees  of  freedom  and  large 
bandwidth.  Isoparametric  modeling  techniques  and  state-of-the- 
art  numerical  solution  methods  are  combined  in  MAGNA  to  provide 
effective  analytical  capabilities  for  finite  strains,  arbitrary 
rotations,  and  elastic-plastic  behavior.  Both  static  and 
transient  dynamic  solution  options  may  be  performed  with  tht? 
program,  as  well  as  natural  frequency/mode  shape  calculations. 
Features  such  as  user  subroutine  interfaces,  post-analys is 


A, 
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graphics,  and  analysis  restart  capabilities  are  included  in 
MAGNA. 

For  the  problem  investigated  herein,  a  two-dimensional 
large  displacement  element,  referred  to  in  MAGNA  as  Element  #9, 
has  been  used  following  a  modification  for  obtaining  plane- 
strain.  This  modification  required  a  user  subroutine  which  is 
presented  in  Appendix  B. 

The  MAGNA  finite  element  program  is  operational  on  the  CDC 
6000  series,  CYBER-74,  and  CYRER-175  computers  with  cupport  CCL 
(CYBER  control  language)  procedures  and  the  segmentation  loader. 
An  example  of  the  job  control  language  (JCL)  used  for  executing 
the  program  during  the  investigations  reported  herein  is  pre¬ 
sented  in  Figure  A-2.  This  JCL  attaches  permanent  files  which 
include  the  input  data  (MAGNADATA) ,  the  user  subroutines 
( TNMAGNA) ,  and  a  dynamic  analysis  restart  file.  It  also  forms  a 
post-processing  file  for  plotting  the  plate  response  (Tape 
L06296)  and  forms  the  next  dynamic  analysis  restart  file  (Tape 
L06091)  following  execution  of  the  program. 

The  nonlinear  dynamic  response  investigation  was  started  by 
first  performing  a  static  displacement  analysis  to  obtain  an 
initial  displacement  (ID)  for  the  initial  transient  dynamic 
analysis.  In  the  time  plots  of  the  section  entitled,  "Analysis 
Results",  time  zero  is  indicated  by  t  =  2.0  because  of  the  static 
analysis  increments  required  for  a  converge  solution.  The  first 
dynamic  analysis  used  the  static  restart  file.  After  a  specified 
number  of  time  increments,  a  dynamic  restart  tape  was  formed  for 
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100=NT  2  j  T1200 , 1 0  3200 ,  CM1O500  ,  GE2  D800956, NOLL, 56832 

110=SET=R1=MFL 

120=ATTACH , A , MAGNAD ATA ,  C  Y =5 

130=ATTACH , B ,TNMAGN A, CY=3 

140=REWIND,A,B 

150=COPY,A,TAPES 

1 60=COP Y  y  B  y  US  RS  UB 

170=RETURN,AyB 

1 80=REW I ND  >  T A  PE S , US  RS  UB 

190=ATTACH , OLDTAP, TDF1 , CY=2 

200=REWIND,OLDTAP 

210=SKIPKyOLDTAPylly0yB 

220<OPYBR ,  0  LDTAP ,  TAPE23 

230=REW I ND ,T APE23 

240=RETURN,OLDTAP 

250=REQUEST,MPOST,GE,R1NG,VSH=L06296 
260=REQUEST>NRSTAP,GE,RING,VSN=L06091 
270=ATTACH,P,MAGNAJCL,ID=BROCKMAN,MR=1 
280=BEGI  N  ,XMAGNA,  P„  USRSUB,  Rl+B 
290=UNLOAD,MPOST,NRSTAP 


Figure  A-2.  Job  Control  Language  for  Executing  MAGNA 


the  next  dynamic  analysis.  The  results  of  each  analysis  were 
placed  on  a  post-processing  file  for  plotting. 

The  dynamic  restart  file  was  then  used  to  continue  the 
analysis,  again  forming  both  post-processing  and  restart  files. 
This  procedure  could  be  repeated  until  the  panel  response  versus 
time  was  constant  (limit  cycle  had  been  reached).  Figure  A-3 
presents  a  procedure  for  using  the  post-processing  data  file  to 
obtain  panel  nodal  point  displacements  plotted  versus  time. 

Figure  A-4  presents  typical  output  results  using  the  FDL/FIBRC 
IMLAC  in  a  Tecktronic  Mode  of  generating  one-line  plots  of  panel 
response . During  this  investigation,  all  MAGNA  computer  runs  and 
plotting  were  accomplished  using  terminal  inputs. 

Piston  Theory  Aerodynamics 

Piston  Theory  aerodynamics  was  used  in  this  study  to  load 
the  panel.  For  Piston  Theory,  the  unsteady  aerodynamic  pressures 
are  defined  as: 


A  ( x ,  t ) 


•  • 
_  -2q  3w  _1  3_w 

8  lax  U  3t_ 


(A-3) 


Using  the  nondimensional  aerodynamic  parameters  defined  in 
References  1  and  2,  that  is: 


\  =  Ifla: 

A  3D 


and  m  fh 

m 


Equation  (A-3)  can  be  expressed  as: 
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100=. PROC, NOLL. 

110=ATTACH,PROCFIL,ID=D800236. 

120=BEGIN,NOSF1LE. 

130=GET, WRTFI L, ID=D800236 . 

140=FTN,I=WRTFIL,L=0. 

150=ATTACH,TAPE99,POSTPR,CY=1. 

160=LGO . 

170=RETURN , PROOF I L> WRTFI L^ LGO , TAPE  99 . 
180=FET,XYPLOT, ID=D800236 . 
190=FTN,I=XYPLOT,L=0. 
200=ATTACH,LIB1,TEKLIB,ID=LIBRARY,SN=ASD. 
210=ATTACH ,LIB2,P L0T3D > ID=KI NG . 
220=LIBRARY,LIB1,LIB2. 

230=LGO . 

240=RETURN,LGOAIB1,LIB2,XYPLOT,TAPE2,TAPE30. 

250=REVERT. 

260=XEOR 

270=XEOF 


Figure  A-3.  Procedure  File  for  Executing 
MAGNA  Plot  Capability 
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BEGIN,  NOLL,  PROCNOLL 
AT  CV=»331  SN=AFFDL 
?FN  IS 
PROCFIL: 

AT  CV»001  SN=AFFDL 

NOSFILE  VERSION  3  READY. 

OTHER  PROCFIL  OPTIONS 
EDTFILE ,  ALLFILE,  TEKFILE,  HPFILE 
FILE  NAME  WRTFIL  HAS  BEEN  RETRIEVED 
.363  CP  SECONDS  COMPILATION  TIME 

MAGNA  POST  FILE  TRANSLATOR  FOR  XYPLOT 

WRITE  DISPLACEMENT  OR  STRESS  DATA  ID, 9) . D 

HOW  MANY  CURVES  WILL  BE  PLOTTED ?1 
ENTER  NODE  NUMBER  AND  DISPLACEMENTSl . 2 

FOR  CURVE  NO.  1  WANT  THE  FIRST  POINT  TO  3E  AT  J.3.3.37N 
STOP 

333330  MAXIMUM  EXECUTION  FL. 

3.371  CP  SECONDS  EXECUTION  TIME. 

FILE  NAME  XYPLOT  HAS  SEEN  RETRIEVED 

PLOT  ON  HP  OR  TEKTRONIX  (l»HP,2=TEK) . : 

UNSATISFIED  EXTERNAL  REF  —  SETIN 
NON- FATAL  LOADER  ERRORS  - 
UNSATISFIED  EXTERNAL  REF  —  SETUU  2 
READ  DATA  FROM  FILE  OR  TERMINAL  < 1 =FILE , 2 =TERM ! . : 1 

NEED  CONVERSION  OPTION  FOR  DATA  SET  NO.  1? . :N 

HOW  MANY  CURVES  DO  YOU  WISH  TO  PLOT . :1 


WHICH  SETS  OF  DATA  DO  YOU  WISH  TO  PLOT.. 
LABELS  FOR  THIS  PLOT  ARE  COMPUTED  TO  SE: 
X-MAXIMUM  =  2.08 

Y -MAXIMUM  *  .34 

X-MINIMUM  =  1.99 

Y -MINIMUM  -  -.05 

CHANGE  THESE  VALUES  (Y,N> . 


ENTER  MAXIMUM  X-LABEL. . 
ENTER  MAXIMUM  Y-LA3EL. . 
ENTER  MINIMUM  X-LABEL.. 
ENTER  MINIMUM  Y-LABEL. . 
INTERVAL  SIZE  ON  X-AXIS 
INTERVAL  SIZE  ON  Y-AXIS 


ENTER  X-AXIS  LABEL,  (MAX.  43  CHAR.) . :  TIME, SEC 

ENTER  Y-AXIS  LABEL,  (MAX.  43  CHAR.) . :  NODAL  ?T  51  DISPLACEMENT ,  IN 


FOR  DATA  SET  NUMBER  1 
LINE  OPTIONS  ARE: 

SOLID  LINE  CURVE,  NO  SYMBOL,  TYPE  -1 
DASH  LINE,  NO  SYMBOL,  TYPE  -2 
SYMBOL  AT  EACH  POINT,  TYPE  1 
SYMBOLS  CONNECTED  WITH  SOLID  LINE,  TYPE  2 
SYMBOLS  CONNECTED  WITH  DASH  LINE,  TYPE  3. 
ENTER  CURVE  TYPE . 


Figure  A-4.  Example  Output  of  MAGNA  Plot  Routine 
Using  Procedure  File 
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A  (  x  ,  t ) 


-XD  /sw  /a  ^1^Mpmh 
a3  l  3x  /  3 DX 


( A-4  ) 


By  grouping  the  dimensional  and  material  properties  of  the  panel, 
the  aerodynamic  pressure  can  be  represented  as  a  function  of 
only  X,  u/M,  and  M  as: 


A  ( x ,  t ) 


C,X 


3w 
3  X 


+  c1c2 


( A-  5 ) 


and 


Using  an  aluminum  plate  (p^  =  — .J-  — ■ 5e-?  and  E  =  107  psi) 
and  assuming  that  the  plate  stiffnes,  _D,  equals  Eh3/12,  the 
aerodynamic  pressure  becomes: 


A  ( x  ,  t ) 


- . 0225X 

3x 


.0001322 


3  w 
3  t 


( A~6 ) 


Now  it  is  necessary  to  define  the  nodal  point 
slopes,  3w/3x,  and  velocities,  3w/3t,  for  Equation  (A-6).  The 
nodal  point  displacements  in  an  element  local  coordinate  system 
(referring  to  the  sketch  below)  are  defined  as: 

w(£)  =  niwa  +  N2WB  +  ^3WC  *  N^W  (A-7) 
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I 


-1 

xe  * 

►-= - « 

o 

4 

element 


xe+1 


The  shape  functions  in  Equation  (A-7)  are  defined  as: 

Nx  =  -U/lMl-O;  N2  =  (1-S2);  and  N3  =  (C/2H1  +  C)  (A-8) 
These  shape  functions  provided  parabolic  displacement  approxima¬ 
tion  across  the  element.  Since: 


-  r VS+i.  +  rf 

x  (  2  J  +  ^  l 


Xe^rle,  .  ~  +  I 

2  J  xe  2 


(A-9) 


then  the  slope  of  a  node  becomes: 


I? |o  *  4^  *  (“a  -  2wb  -  ”c>4 


and  the  node  velocity  becomes: 


Tt  U)  =  NlWA  +  N2WB  +  N3WC 


( A-10 ) 


( a-ii: 


From  the  Principle  of  Virtual  Workf  the  work  done  by  the  aero¬ 
dynamic  pressure  A(x,t)  was  equated  to  the  work  done  by  the  nodal 


■?  \ 


forces  moving  through  the  nodal  displacements.  In  other  words, 
for  a  particular  element,  the  work  done  is: 


X  ,  , 

e+1 

f  \ 

fa 

T 

/  \ 

WA 

W  =  /  A(x,t)w(x)dx= 

< 

fb 

WB 

X 

e 

iFcJ 

,wc  . 

or 


1 

f  1 

fa 

T 

'  WA  ) 

W  =  V2  /  A(?,t)w(c  )d£  =  < 

fb 

»  < 

WB 

-1 

,Fc. 

,wc  J 

K1 

I  Pb  i 

l  pc  J 

and 

PA  =  A( — 1 , t ) , 

PB  =  A( 0 , t ) , 

Pc  =  A( 1 , t )  for  the  eth  element, 
substitution  into  Equation  (A-12)  gives: 


With  w(^)  =  N 


w. 


w. 


w. 


A(?  ,t)  =  NJ 


(A-12) 


N 


or 


w 

w 

w 


A 

B 

C 
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{ A-l 3 ) 


fbi-Mp 

V  F_ 


/  NNTdC 


-1 


By  evaluating  the  integral  for  the  shape  functions  N-^,  1^2 ,  and 
N-j,  Equation  (A-13)  becomes: 


v 


B 


’c 


1 

60 


r  8 

4 

u  -2 


4 

32 

4 


-2  “I 

4 

8  -* 


V 


( A- 14 ) 


Now  using  the  nodal  point  displacement  and  velocities,  the 
pressures  evaluated  at  the  element  nodal  points  become: 


PA  =  C1X 


<-3wa  +  4wq 


wc)  +  C2J  w; 


( A-15a ) 


PB  =  C1A 


Wc  -  WA  +  C2/—W  WF 


( A-l 5b ) 


pc  -  cix 


(WA  -  4wB  +  3V  +  C2v 


. 

3X  WC 


(A- 15c) 


By  substituting  these  expressions  into  Equation  (A-14)  we  obtain 
an  equation  which  now  relates  the  nodal  point  forces  to  the  nodal 
point  displacements  and  velocities.  Equation  (A-16)  becomes: 
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-  4 

2 

-1  - 

U 

A 

2 

16 

2 

' 

-  -1 

2 

4  - 

.  w~  - 

(A-16) 


At  common  nodal  points  between  two  elements,  the  nodal  forces  are 
added.  Equation  (A-16)  was  programmed  in  the  ULOAD  subroutine 
for  calculating  the  nodal  point  forces  based  on  Piston  Theory 
aerodynamics.  Appendix  C  presents  a  listing  of  ULOAD. 

Analysis  Results 

This  section  presents  the  results  of  a  vibration  analysis 
and  a  nonlinear  static  analysis,  and  briefly  summarizes  the 
results  of  the  nonlinear  transient  dynamic  analyses.  Some  of  the 
response  data  provided  herein  are  sketches  of  (w/h)  versus  a  non- 
dimensional  time,  t,  defined  as: 

t  =  t(D/pmha4)/2 

For  the  plate  used  in  this  study,  t  =  17.02t.  These  sketches 
were  made  prior  to  obtaining  on-line  plot  capability.  The 
sketches  are  somewhat  rough  in  that  there  was  no  attempt  to 
obtain  an  accurate  time  history;  only  the  peaks  and  zeros  of  the 
time  history  were  plotted.  All  on-line  plots  present  _w  versus 
real  time,  t_,  with  time  zero  beginning  at  t  =  2.0  secs  as 
described  earlier.  The  time  increment  jump  from  0.0  to  2.0  was 


41 


caused  during  the  generation  of  the  static  restart  tape  which 
required  two  increments  for  convergence  (2.0  secs).  Therefore, 
the  first  dynamics  run,  using  a  static  restart  file,  would  begin 
at  t  =  2.0  secs. 

Vibration  Analysis 

The  eigenvalue  solution  option  of  MAGNA  was  used  to  obtain 
the  first  five  mode  shapes  and  natural  frequencies  of  the  plate. 
These  data  are  found  in  Table  A-l.  The  mode  shapes  are  presented 
for  the  nodal  points  through  the  center  of  the  plate  (nodes  #45- 
#54  and  the  simply  supported  edges  #1  and  #44).  A  consistent 
mass  representation  was  selected  for  this  analysis.  MAGNA  uses  a 
vector  iteration  procedure  for  obtaining  the  desired  eigenvalues 
and  eigenvectors. 

For  a  simply  supported  beam,  the  natural  frequencies  of  the 
beam  can  be  represented  as: 

a  =  (  nir/a)  2(Eh2/12pm)1/2  ,  n  =  1,2, _ 

This  expression  gives  26.73  Hz,  106.92  Hz,  and  240.57  Hz  for  the 
first  three  natural  frequencies.  These  results  are  in  good 
agreement  with  the  ttequencies  found  in  Table  A-l.  All  the  beam 
results  are  lower,  as  expected.  The  mode  shapes,  which  are  pre¬ 
sented  with  only  three  decimal  digits,  all  appear  reasonable  with 
respect  to  the  mode  shapes  of  a  simply  supported  beam. 

Nonlinear  Static  Analysis 

The  static  analysis  input  data  to  MAGNA  are  presented  in 
Figure  A-5.  The  static  analysis  was  required  to  obtain  an 
initial  displacement  (ID)  for  the  nonlinear  transient  dynamic 
analyses.  Four  initial  displacements  were  investigated  during 


TABLE  A- 1 

NATURAL  FREQUENCIES  AND  MODE  SHAPES 


Mode  No. : 

1 

2 

3 

4 

5 

Freq  (Hz) 

28.04 

112.13 

252.69 

451.18 

711.25 

Modal 

Mode  Shapes 

PT 

1 

2 

3 

4 

5 

1 

.000 

.000 

.000 

.000 

.000 

45 

-.309 

.588 

-.809 

-1.000 

-1.000 

46 

-.588 

.951 

-.951 

-.618 

.000 

47 

-.809 

.  951 

-.309 

.618 

1.000 

48 

-.951 

.588 

.588 

1.000 

.000 

49 

-1.000 

.000 

1.000 

.000 

-1.000 

50 

-.951 

-.588 

.588 

-1.000 

.  000 

51 

-.809 

-.951 

-.309 

-.618 

1.000 

54 

-.  707 

-1.000 

-.706 

.000 

.698 

52 

-.588 

-.951 

-.951 

.618 

.000 

53 

-.309 

-.587 

-.809 

-1.000 
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Figure  A-5.  MAGNA  Data  Input  for  Nonlinear  Static 
Load  Analysis 
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the  course  of  this  study.  The  static  load  and  the  resulting 
plate  displacements  that  were  used  are  shown  in  Table  A-2.  For 
these  analyses,  a  static  load  was  placed  at  the  center  of  the 
plate  (nodal  point  #23).  A  nonlinear  large  displacement 
formulation  was  assumed  in  these  analyses.  For  the  results 
defined  in  the  section  entitled,  "Nonlinear  Transient  Dynamic 
Analyses",  all  time  histories  have  a  specific  initial  displace¬ 
ment  associated  with  it.  These  ID's  are  defined  in  Table  A-2. 

For  a  simply  supported  beam,  the  deflection  at  the  center  can  be 
defined  as: 

w 49  =  Fa3/4Eh3 . 

This  expression  gives  results  which  are  in  good  agreement  with  ID 
#1  of  Table  A-2.  The  other  loads  result  in  MAGNA  predicted 
deflections  which  are  considered  to  be  "large  deflections".  The 
beam  expression  was  derived  assuming  small  displacements,  there¬ 
fore,  correlation  with  ID  #2,  ID  #3,  and  ID  #4  would  be  expected 
to  worsen. 

Nonlinear  Transient  Dynamic  Analyses 

The  nonlinear  transient  dynamic  analyses  were  conducted  for 
three  different  Mach  numbers;  3.0,  1.2,  and  0.8.  When  the  Mach 
number  of  0.8  was  used,  8  was  redfined  as  (1-M^)^2.  Figure  A-6 
presents  a  sample  of  the  data  input  to  MAGNA  for  these  analyses. 
Any  variations  to  X,  y/M,  or  Mach  number  were  accomplished  in  the 
ULOAD  suoroutine. 


TABLE  A- 2 

INITIAL  DISPLACEMENTS  USED  FOR 
TRANSIENT  DYNAMIC  ANALYSES 


Nodal 

PT 

.00357 

Static  Load  at  PT  23  (lbs 

.00714  .01071 

) 

.02570 

ID#1  (in) 

.  ... 

ID#2  (in)  ! 

ID#3  (in) 

ID#4  (in) 

1 

.  0 

.0 

.0 

.0 

45 

-.00087 

-.00162 

-.00224 

-.00396 

46 

-.00166 

-.00311 

-.00430 

-.00762 

47 

-.00232 

-.00434 

-.00601  i 

-.01070 

48 

-.00276 

-.00517  | 

-.00718 

-.01284 

49 

-.00292 

-.00548 

-.00761 

-.01365 

50 

-.00276 

-.00517 

-.00718 

-.01283 

51 

-.00231 

-.00433 

-.00600 

-.01069 

54 

-.00201 

-.00376 

-.00521 

-.00926 

52 

-.00166 

-.00311 

-.00431 

-.00764 

53 

-.00087 

-.00162 

-.00224 

-.00396 
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Figure  A-6.  MAGNA  Data  Input  for  Nolinear  Transient 
Dynamic  Analysis 
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Figures  A-7,  A-8,  and  A-9  present  time  history  plots  of 
nodal  point  #54  (3/4  span)  for  initial  displacements  ID  #1,  ID  #2 
and  ID  #3  and  for  X,  y/M  and  M  equal  to  400,  .01  and  3.0,  respec¬ 

tively.  The  figures  present  nondimensional  displacement,  (w/h), 
versus  nondimensional  time,  t.  As  described  earlier,  the  plots 
are  approximate  and  only  show  trends.  As  can  be  seen  from  the 
figures,  the  displacements  are  approaching  the  thickness  of  the 
plate.  The  frequency  of  the  oscillation  is  about  83  Hz.  Since 
the  objective  of  this  study  is  ultimately  to  investigate  tran¬ 
sonic  effects,  the  studies  at  M  3.0  were  not  continued  to  limit 
cycle. 

Mach  Number  =  1.2 

For  a  Mach  number  of  1.2,  the  aerodynamic  parameters,  X 
and  p/M,  were  selected  to  be  300  and  0.013,  respectively.  Two 
different  initial  displacements  wer  investigated  at  these  condi¬ 
tions.  Figure  A-10  presents  the  time  history  results  for  ID  #4, 
and  Figure  A-ll  presents  results  for  ID  #3.  These  figures  pre¬ 
sent  the  time  histories  of  the  nodal  point  displacement,  _w, 
versus  real  time,  _t_,  for  the  points  through  the  center  of  the 
plate  (points  #45-#54).  Nodal  point  #54  represents  the  3/4  span 
location  of  the  plate  and  is  positioned  correctly  for  increas¬ 
ing  within  the  figures  (between  points  #51  and  #52).  The  nodal 
point  is  described  in  the  ordinate  label  of  each  time  history 
plot.  Also  note  that  in  Figure  A-10,  there  is  a  change  in 
ordinate  scale  after  t  =  2.7  secs. 
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Figure  A-9.  Transient  Dynamic  Response  to  Initial 
Displacement  #3,  PT.54,  (X  =  400, 
y/M  =  .01,  M  =  3.0) 
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After  three  dynamic  restart  cycles,  the  time  history  for  ID 
#4  (Figure  A-10)  reaches  a  limit  cycle  near  t  =  2.21  secs  (.21 
secs).  The  magnitude  of  the  displacement  at  point  #54  is  ±.0144 
inches  or,  in  nond imens ional  form,  (w/h)  =  ±.480.  For  ID  #3 
(Figure  A-ll),  the  time  history  shows  limit  cycle  after  two 
dynamic  restart  cycles.  At  point  #54,  the  displacement  of  the 
limit  cycle  is  ±.012  inches  (w/h  ±.400)  at  t  =  2.13  (.13  secs). 
This  information  indicates  that  the  limit  cycle  amplitude  is 
dependent  not  only  on  Mach  number,  y/M  and  X,  but  also  on  initial 
displacement.  The  maximum  amplitude  of  the  motion  occurs  at 
about  C  =  .7  (70%  span);  the  displacement  then  falls  off  to  zero 
at  nodal  points  #1  and  #44.  The  frequency  of  the  oscillation  is 
about  65  Hz. 

Mach  Number  =  0.8 

Figure  A-12  presents  the  transient  dynamic  analysis  results 
for  a  Mach  number  of  0.8  and  initial  displacement  ID  #3,  again 
with  X  =  300,  and  y/M  =  .013.  The  results  are  similar  to  that 
discussed  in  the  section  entitled,  "Mach  Number  =  1.2".  Motion 
that  was  very  near  the  limit  cycle  was  reached  at  t  =  2.21  (.21 
secs).  The  amplitude  of  the  motion  at  point  #54  was  ±.017  inches 
(w/h  ±.580).  Again,  the  frequency  of  the  oscillation  is  about 
65  Hz. 

A  summary  plot  which  presents  the  deflection  of  the  plate 
during  limit  cycle  for  Mach  numbers  of  0.8  and  1.2,  is  shown  in 
Figure  A-13.  The  response  at  M  =  0.8  is  larger  across  the  plate 
than  the  response  at  M  *  1.2. 


Figure  A-10.  Transient  Dynamic  Response  to  ID#4 
(X  =  300,  u/M  =  .013,  M  =  1.2) 
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Figure  A-12.  (continued) 
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ELEMENT  #9  ROUTINE  MODIFIED  FOR  PLANE  STRAIN 
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APPENDIX  C 

USER  ROUTINE  FOR  PISTON  THEORY  AERODYNAMICS 


1  SUJ:iOUT:f.r  UIOAOCP.WK  ,  TI*'E,OT,NTNODV,NRDCF,NRCCOE,t'lTblS, 

_»NRU,N®V.NTt',,NOUTl 
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j,  -  N&  _  i iSLf.ii* s*4 


1**  _ N3  =  (IELF-i>»«*A 
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1Z LC  =  W< inti 

2*  o<NAt  =  ,  27S*PACTl*<-’.»3I3PA*».*DISP3-0iSPC)* 
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APPENDIX  D 


COMMENTS  ON  NUMERICAL  INSTABILITIES 
WHEN  BOTH  PANEL  AND  AERODYNAMIC  EQUATIONS  ARE  NONLINEAR 

An  analysis  was  initiated  to  study  the  panel  response 
represented  by  a  nonlinear  structural  equation  (Von  Kantian' s 
large  deflection)  and  a  nonlinear  aerodynamic  equation  repre¬ 
sented  by  NLR-LTRAN2  numerical  code.  A  time  dependent  response 
study  requires  the  simultaneous  time  integration  of  the  struc¬ 
tural  and  aerodynamic  nonlinear  equations.  Here,  the  panel 
response  was  represented  by  a  linear  superposition  of  assumed 
modal  terms  even  though  both  the  structural  equations  and  aero¬ 
dynamic  equations  are  nonlinear.  Using  the  Galerkin  procedure, 
the  spatial  dependency  in  the  structural  equations  was  eliminated 
and  were  integrated  simultaneously  with  the  aerodynamic  equations 
represented  by  the  NLR  version  of  LTRAN2 .  The  NLR  version  of 
LTRAN2  included  a  modification  to  account  for  the  changing 
induced  angle  of  attack  resulting  from  the  panel  deformation 
which  was  incorporated  into  a  chord  deformation  option  contained 
internal  to  NLR-LTRAN2. 

The  simultaneous  time  integration  procedure  allows  one  to 
determine  panel  deflection,  velocities  and  accelerations.  Once 
the  initial  conditions  for  the  panel  deformation  are  provided, 
the  effective  induced  angle  of  attack  and  its  time  derivatives 
can  be  calculated  for  each  point  on  the  panel.  The  resulting 
aerodynamic  pressures  can  then  be  calculated.  The  panel  response 
or  deformation  at  the  next  time  increment  can  be  determined  and 


the  process  repeated  for  each  desired  ti.<e  step.  Prior  to 
releasing  the  panel  for  response  to  the  aerodynamic  pressure,  the 
aerodynamic  equations  were  integrated  in  time  for  an  assumed 
simple  harmonic  panel  deformation  until  the  transient  prediction 
of  the  aerodynamic  pressures  disappeared  and  the  pressures  became 
periodic.  This  assumed  simple  harmonic  motion  generally  required 
the  time  for  two  periods  of  panel  oscillation  to  insure  that  the 
convergence  criteria  specified  in  NLR-LTRAN2  was  satisfied. 

For  the  panel  and  airstream  parameters  investigated,  the 
simultaneous  integration  of  nonlinear  structural  and  aerodynamic 
equations  were  unsuccessful  in  obtaining  stable  panel  responses. 
In  all  cases  the  panel  response  diverged  sometimes  quite  quickly 
when  the  panel  was  freed  to  response  to  the  transonic  aerodynamic 
pressure.  The  diverging  panel  response  was  attributed  to  numeri¬ 
cal  instabilities  since  previous  investigations  using  a  linear 
aerodynamic  representation  at  greater  Mach  numbers  (piston 
theory)  yields  a  stable  panel  response  for  certain  selected  panel 
parameters.  Various  magnitudes  of  initial  displacement,  velocity 
and  acceleration,  both  singly  and  in  combination,  were  investi¬ 
gated.  The  time  step  for  the  time  integration  was  varied  to  very 
small  values  within  computational  considerations.  Also,  the 
number  of  cycles  for  forcing  the  plate  motion  in  a  sinusoidal 
fashion  was  varied  to  assure  all  transients  were  small.  In  all 
cases  considered,  the  panel  response  was  divergent  after  a  given 
period  of  time. 
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The  following  suggestions  are  offered  as  possible  courses 
of  the  aforementioned  numerical  instabilities  for  consideration 
by  future  investigators. 

1.  The  chord  deformation  option  within  NLR-LTRAN2  was  a 
constant  source  of  convergence  problems.  Even  for  very  small 
panel  deformations  the  local  angles  of  attack  (slopes  of  panel 
deformation)  can  become  large  enough  to  exceed  the  convergence 
criteria  specified  in  NLR-LTRAN2 .  This  would  result  in  an  over 
prediction  of  aerodynamic  pressure  with  no  warning  from  t!~^ 
computational  code.  This  investigator  is  unaware  of  any  verifi¬ 
cation  that  the  chord  deformation  option  has  been  accomplished 
and  any  coding  errors  would  surely  yield  a  source  of  numerical 
instabilities.  Further  investigation  of  the  aerodynamic 
pressures  predicted  by  this  option  of  NLR-LTRAN2  for  a  given 
chord  deformation  would  seem  to  be  warranted  to  insure  numerical 
stability. 

2.  Representation  of  the  nonlinear  panel  response  to  the 
nonlinear  aerodynamic  pressures  by  a  linear  superposition  of 
natural  modes  may  be  another  source  of  numerical  instability  and 
is  a  somewhat  questionable  approach.  As  suggested  in  Appendix  A, 
the  representation  of  the  large  deflection  structural  equations 
with  a  large  deflection  finite  element  has  some  merit  for  inves¬ 
tigation  of  nonlinear  response. 

3.  In  performing  the  simultaneous  integration  of  struc¬ 
tural  and  aerodynamic  equations  there  is  by  necessity  a  differ¬ 
ence  in  time  step  increments  required  primarily  to  insure 
convergence  of  the  chord  deformation  option  used  in  NLR-LTRAN2. 


The  difference  in  the  time  scaling  between  the  structure  and 
aerodynamic  equations  were  not  fully  investigated  and  would  be 
another  source  of  numerical  instability. 

4.  The  stable  response  obtained  using  a  linear  aerody¬ 
namic  theory  were  conducted  at  supersonic  Mach  numbers  whereas 
the  Mach  numbers  selected  in  this  investigation  were  in  the  so 
called  "sub-transonic"  Mach  number  regime.  The  flow  field 
characteristics  are  radically  different  than  the  supersonic  Mach 
regime  and  would  require  numerical  carefulness  to  insure  numeri¬ 
cal  stability.  Verification  of  the  possibility  of  stable 
responses  in  this  "sub-transonic"  region  needs  to  be  established. 

5.  Another  possible  cause  could  be  the  time  step  size 
resulting  in  an  instability  generated  by  the  motion  of  a  shock 
wave  due  to  the  mixed  numerical  differencing  in  NLR-LTRAN2 .  It 
would  be  necessary  to  choose  the  time  interval  small  enough  such 
that  shock  waves  do  not  travel  more  than  one  mesh  point  in  the  x 
direction  over  a  single  time  step.  However,  very  small  time 
steps  may  result  in  prohibitive  computational  times. 

6.  Finally,  the  initial  panel  deformation  and  aerodynamic 
pressure  generated  initially  before  release  of  the  panel  to 
response  to  the  aerodynamic  pressure  needs  a  more  complete  inves¬ 
tigation  when  both  nonlinear  structural  and  nonlinear  aerodynamic 
equations  are  integrated  simultaneously. 


